#pragma once
#include <memory>
#include <string>
#include "RI.h"
#include "Grid.h"
#include "MVoperator.h"

using namespace std;
class Cycle
{
private:
    // 执行一次Jacobi加权迭代
    void Jacobi(vector<double>& u, vector<double>& f) {
        int n = u.size() + 1;
	vector<double> cf = gM.M_times_V(gM.get_D(n),f);
	for (int i = 0 ; i < cf.size(); i++)
	{
	    cf[i] = (double)2*cf[i]/3;
	}
        u = gM.V_plus_V(gM.M_times_V(gM.get_Tw(n),u),cf);// u = Tw * u + cf
    };
    // 计算残差
    void compute_residual(vector<double>& u, vector<double>& f, vector<double>& r) {
        int n = u.size()+1;
        r = gM.V_minus_V(f, gM.M_times_V(gM.get_A(n), u));// r = f - Au
    };
    InterpolationOperator* I;
    RestrictionOperator* R;
public:
    Cycle(InterpolationOperator* I, RestrictionOperator* R): I(I), R(R){};
    MatrixVectorOperator gM;
    vector<double> Vcycle(vector<double> &v, vector<double> F) // 进行V-cycle
    {
	int N = F.size() + 1;
	//cout << v.size() << " " << F.size() << endl;
	Jacobi(v, F);
	//cout << N << endl;
	vector<double> r(F.size()),e(F.size());
	if (N != 2)
	{
	    compute_residual(v, F, r);
	    vector<double> r_coarse(F.size()/2), e_coarse(F.size()/2);
	    r_coarse = R->restrict(r);
	    Vcycle(e_coarse, r_coarse);
	    e = I->interpolate(e_coarse);
	    v = gM.V_plus_V(v,e);
	}
	Jacobi(v, F);
	return v;
    };
    vector<double> FMG(vector<double> &v, vector<double> F) // 进行FMG
    {
	int N = F.size() + 1;
	// vector<double> v(N-1);
	if (N != 2)
	{
	    vector<double> F_2,v_2(N/2 - 1);
	    F_2 = R->restrict(F);
	    v_2 = FMG(v_2, F_2);
	    v = I->interpolate(v_2);
	}
	Vcycle(v, F);
	return v;
    };
};
